A procedure to produce excess, probability and significance maps and 
to compute point-sources flux upper limits 

P. Billoir (Billoir@in2p3.fr) & A. Letessier-Selvon (Antoine.Letessier-Selvon@in2p3.fr) 

A short note to propose a procedure to construct excess maps, probability maps and to calculate point source 
flux upper limits. 



1. Excess maps 

In the following we do not discuss how the coverage map is constructed, be it with shuffling, with the semi- 
analytical method or whatever we simply assume that we have, in the case of a perfectly isotropic CR sky, 
a prediction for the average number of cosmic ray we should observe from a given direction on the sky. Let 
Mbg{fl) be this coverage function and Mbg{k) its value integrated over the pixel k (centered on direction flk) 
of a pixelized representation of the sky. 

Given a set of events and their individual pointing accuracy we can construct a CR density map AI {k) just by 
counting the number of events that fall in pixel k. We can also construct a smooth CR density map distributing 
those events over the map according to their pointing accuracy. Let Alg be such a smooth function and Ms(k) 
its value integrated over pixel k. 

We can construct excess map or relative excess map as : 

1. Me(fc) = Af,(fc)-Mb<,(fc);or 

2. Mreik) = Ms{k)/Mbcj{k). 

Since we cannot predict from the exposure (we do not know the flux) the total number of events expected in 
case of an isotropic sky we normalize Mi,g to the total number of events observed {N = J Ms{fl)dil). It 
can be useful to filter (convolute) the above maps with a global Gaussian point spread function of the type 
P{Q,9Q)dfl = /(^^o)rfJ7 where 9 is the angular distance to the center of the pixel we are filtering, to 
emphasize (or to enhance) structures that may be present at a given scale ^o- 



2. Probability maps 



In a probability map, each pixel carries the probability that a uniform sky leads to a larger number of event 
observed in that pixel than what we have measured. Without any filtering and if we do not distribute events on 
the sky according to their pointing accuracy, constructing such a map is straightforward as the above probability 
is given by Poisson statistics. 

Mprobik) = Pin > M{k)-Mba{k)) (1) 



If we have done filtering we can either use a Monte Carlo technique to estimate the above probability or 
approximate the filtered random variable distribution with an appropriate normal distribution whose parameters 
are derived below in the case of a Gaussian filtering. 
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Figure 1. True background distribution after Gaussian filtering on a 1.5° scale (size of filtering does not matter up to at 
least 15°). The average number of events (Neo) in a 1.5° disk is 5 in the left plot and 50 in the right one. The overlayed 
function is the approximated normal distribution {J\f(2Neo, Neo)) with parameters as predicted by Eq. 3 and Eq. 4. 

2.1 MC calculation 



For the MC technique we can construct smoothed density map M]""^ in exactly the same way as we have 
constructed our signal density map Ms- The corresponding probability map will then be given by 



Mprob{k) 



number of maps satisfying: {M™'^(fc) > Ms{k)} 
N 



(2) 



2.2 Gauss approximation 

For Gaussian filtering with parameter Oq << 1, such that we can work on the plane in a disk of radius 6 
ignoring the contribution of the points outside of this disk and the variation of the coverage on that disk, we 
can model the distribution of the expected background. Let /i — N/{ttQ'^) be the uniform CR background 
density then the average weight < w > of a background event is : 

<w>= ^27r£ 0e-''/(''o)d0 = 2 (I)' = 2so (3) 

The above integral is still accurate to a few percent for 6q as large as 15° and Q — 60°. 
While the variance of the weights is given by : 

Viw] =<n? > - <w eo(l - 4£o) ^ eo (4) 



since we assume that eq << 1- 

For N events the random variable X=^^ Wi has an average <X> =2Neo and a variance V[X] = Neq. The 
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probability density function of X can be adequately modeled (especially on the excess side) by a normal law 
J\f{2Neo , Neo ) even for an average background count in a circle of radius eq {Neo) as low as 5 . See Figure 1 . 



The probabiUty map is then obtained from the filtered signal and coverage map as : 

'Ms{k) 

Or equivalently a significance map can be obtained as : 



r+oo 

Mprob{k)= Af{x;M^'{k),M^%k)/2)dx (5) 

JMJk) 



M.{k)-Mrik) 
V Mr-(fc)/2 

Remark : If the number of pixels in the maps are very large so that the probabiUty to have an event falling in 
a given pixel is much smaller than 1 then, without filtering, the probability map gets very difficult to interpret 
(at least visually) as it wiU be only composed of pixel with values or 1. 



3. Optimizing point source detection 

Here we want to address the question of how to optimize the search region to detect an eventual point source 
located at ilo. For simplicity we rotate our coordinate system so that f^o = 0. In the following we assume that 
the aperture of the array does not vary significantly on the scale of the detector angular resolution (a), therefore 
the background can be considered uniform in the region we are looking at. 

The average number of events from a uniform CR distribution expected from that direction is given by : 

^lbg = I A{e)^cr{e) I w{n)dndE (?) 

JAE J An 

where ^cr{E) is the uniform CR flux per unit solid angle , time, surface and energy. A{E) is the aperture 
(per unit time, surface, and energy) which we assumed independent of O, and 1^(0) is the weight function we 
want to optimize to detect point-like sources. 

From an hypothetical point source s in that same direction the average number of events is given by : 

IJ.s= f A{E)^,{E) [ W{fl)r]{fl)dmE (8) 
Jae Jaq 

where ^s{E) is now the source flux per unit time surface and energy and ri{Cl) is the detector point spread 
function (PSF), assumed independent of energy for simplicity. 

The significance S of the signal is given by the ratio of the signal from the source to the background 
fluctuation i^yv[nbg\). To maximize S we need to maximize this ratio at each energy that is to maximize : 

R{W) = -^ f" ^ (9) 

^Jv [J^^W{n)dil] 

with respect to W{il). 

A general result of signal processing theory tells us that the optimal choice for W is the detector response 
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Figure 2. Ratio of Roausa to R'p (solid-line), Rcauss (dotted-line) and V[ns]/<ns> as a function of a. 

function, i.e. taking W proportional to ?7(f2) in our case. With such a choice and using W{^) = e"^^/^^'^^^ 
the i?-ratio becomes : 



Rca 



(10) 



where we have used Eq. 4 for the variance of the background. 

It has been argued that a top-hat weight function W can also maximize this ratio. In case of a top-hat weighting 
function Wa {^) = H(0 — Pa) where H is the Heavyside function and /? the distance in units of a up to which 
we want to extend the source integration, we have : 



^0 = / „^ = ^ I a I ^Ga 



^/n(3^a^^CR{E) 



(11) 



Rp reaches a maximum value of 0.9 x Rcauss for = e^^/^ — 1 (i.e. j3 = 1.585). Which shows that a top-hat 
weighting function is always less efficient than weighting with the detector PSF. 

We studied the evolution of the R-ratio in the case where the PSF used for filtering has a different parameter 
than the true PSF of the detector. Calling a the ratio between the filter parameter and the true PSF parameter 
Eq. 10 and Eq, 1 1 become : 



R, 



2a 



Gauss 



1 + a 



_ d1 

2 -^Gauss 



/ 1 _ e-("/3)V2 \ ^ 
and R% = 2\ 1 i?^ 



Gauss 



(12) 



We show on Figure 3 the evolution of Rcauss ^ function of alpha as well the ratio of Roauss ^0 
P = 1.585(7. We see that for all a Rcauss l^rg^r than i?^ and, as expected, that -Rcauss maximum for 
a = 1.0. 
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4. Flux upper limits 

To obtain flux upper limits at a certain confidence level (CL) one needs first to derive from the number of events 
observed and the expected background contribution an upper hmit, at the same CL, on the source contribution 
to the observation. 

4.1 Top-Hat weighting functions 

If we use a top-hat weighting function to count our events around a certain direction then the source contribution 
upper hmit at CL a, in that window can be directly derived using Poisson statistics. We call /x" this upper limit 
and we have': 

where ni,g is the expected number of count in our top-hat window and iJ.a = /J^" + ntg is such that the Poisson 
cumulative probability of observing more than Hobs events in the window is larger than a : 

P{n > Hobs; IJ'a) = X] P{i;iJ'a) > a. (13) 

or, equivalently such that the probability to observe at most Uobs events is less than 1 — a: 

nobs 

P{n<=nobs;iJ'a) = ^P{i;i^a) <'i--0!. (14) 

This formula can be solved numerically or can be approximated for iibg sufficiently large (above about 10 for 
CL in the range [0.8-0.95], larger CL would require larger average background) using a Gaussian distribution 
iJ\f{nbg, Ubg) for the background contribution. This leads to : 

= (2no6, + Cl + ^^UobsCl + Ci)l2 - Ubg (15) 

where C„ is such that /f^ 7V(a;; 0, l)dx = a (C„ = 1.64 for a=95%). 

4.2 Gauss weighting functions 

In this case the event count around a given direction does not follow a Poisson distribution (it is no longer 
integer) so we can either rely on a Monte Carlo evaluation of /i", throwing background events uniformly 
around the source direction with an average background count equal to our background expectation in the 
filter window and throwing source events with a distribution proportional to / ^^'^ ' finding the source 
normalization for which in CL% of the case the MC count is larger than the observed count^. 

Alternatively we can model the distribution of the sum of the background and signal count and derive the upper 

^ Formally the variable /i™ , n[,g and rifjg should all carry the letters TH to mark the fact they are all computed using a Top-Hat window 
filter. However, For simplicity we dropped this symbol. 

^Here also the variables should carry the letter G to remind they are now computed with a Gaussian filter but again we have dropped 
this superscipt for clarity. 
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limit /z" from this distribution. For the background we have already shown that the count distribution can be 
approximated by a normal law Af{nbg , nbg/2) a similar calculation gives for the signal weights : 

<w>= r^27r / e-^'/^''''^0de = ^ (16) 



and 



< yj^ >= J—27T I e-3^'/(2-')^rf^ = - (17) 
ZTTCT^ J 3 

For a source of average intensity Vs our event count is rig = X^^Li where n is distributed according to a 
Poisson law of average Ug so our average count becomes : 

n 

^P(n;z^s)^< Wi>=y (18) 

n i=l 

with variance : 

n n „ 

V[ns] = V{n; u,) < ^ u;. ^ > -(i/, < to >)2 = ^ = ^ (19) 



Using for the distribution of the total count a Gaussian approximation jV(nfeg + fis,nhg/2 + 2/is/3) the upper 
limit of CL a on the source count is given by the solution of 



riobs - (ribg + IJ's) ^ Ca\J^ + ^ (20) 
which can easily be solved analytically. 

On Figure 3 we show the evolution of the ratio V[ns]/ns as a function of the a parameter introduced in 
section 3 V[ns]/ ij,g{a) = (1 + Q;^)/(l + 2a^). We see that for a > 1 Eq. 19 is an upper bound on V[ng] and 
will therefore lead to conservative upper limit if we mistakenly use a filter with a larger a parameter than the 
detector PSF 

In Table 1 we compare the 95% CL upper limit obtained from Eq. 13, 15 and 20 for various conditions of 
background and for an observed count equal to the background prediction. For the approximation of Eq. 20 we 
indicate the true CL of this upper limit as extracted from a Monte Carlo simulation. 



4.3 Fluxes 

If the source spectral shape is the same as the overall CR spectral shape in the domain AE where we want to 
place a flux upper limit, that is if $s(-E) oc ^cr{E). then the aperture part of the integral contribute in the 
same way for the background and for the signal and we can relate the source flux upper limit to the ratio of /u" 
to our expected background in the f oUowing way : 

^ 1-6-/5^2 ntgiTH) mg{TH) ^ ' 

for an optimal top-hat window of width /3 = 1.585 in units of a the detector Gaussian PSF parameter, and, for 
a Gauss filter of parameter a. 



Excess maps, source flux upper limits 1 



ribg (in a l-cr disk) 


2.5 


5.0 


10.0 


25.0 


50.0 


Top-Hat window (exact) 


5.6 


8.1 


9.8 


14.9 


20.5 


Top-Hat window (approx.) 


5.3 (0.94) 


7.8 (0.94) 


9.5 (0.95) 


14.6 (0.95) 


20.2 (0.95) 


Gauss filter (exact) 
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Gauss filter (approx.) 


3.6 (0.96) 


4.7 (0.96) 
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9.1 (0.96) 
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Table 1. (95% CL upper limit on the source count for an optimal top-hat window and a Gaussian filter with the detector 
PSF parameter. In all case the observed number of events was taken to be equal to the background expectation (or the 
nearest interger if needed). The exact line corresponds to the exact solution given by the numerical solution or an MC 
integration. The approx. line corresponds to the Gauss approximation with its corresponding CL (as computed from a MC) 
in parenthesis, 

where we have explicitely indicated that the values of the source count upper-limit and of the expected back- 
ground count depend on the filter type used. 

Given the values presented in Table 1 , the upper limits from the two methods only vary by a factor 

/ = (4/32^^(G)/(3.5/.:(ri?) (23) 

which is approximately 15% giving a small advantage to the Gauss filter. We should however keep in mind 
that this is the best case comparison for the top-hat window (1.585c7) as any other value (possible if one does 
not know exactly the detector resolution) would increase the difference. 

5. Conclusions 

We have shown that using a Gauss filter to produce sky maps or to compute flux upper limits is an optimal 
choice to enhance the sky feature at a particular scale. However, the gain in sensitivity is rather small (15%) 
compared to a top-hat filter, if the detector resolution is properly know. In other cases the improvement offered 
by Gauss filtering is always larger. For situation where the angular resolution varies from event to event, we 
have shown that using the largest value give good results and safe upper limits. It is however better, and all 
the above considerations are still valid, to adapt the filter on an event by event basis. That is, for each event 
entering the count, to filter it with a Gauss filter matching its resolution. 



